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A detailed description is provided of a new Worm Algorithm, enabling the accurate computation 
of thermodynamic properties of quantum many-body systems in continuous space, at finite temper- 
ature. The algorithm is formulated within the general Path Integral Monte Carlo (PIMC) scheme, 
but also allows one to perform quantum simulations in the grand canonical ensemble, as well as to 
compute off-diagonal imaginary-time correlation functions, such as the Matsubara Green function, 
simultaneously with diagonal observables. Another important innovation consists of the expansion 
of the attractive part of the pairwise potential energy into elementary (diagrammatic) contributions, 
which are then statistically sampled. This affords a complete microscopic account of the long-range 
part of the potential energy, while keeping the computational complexity of all updates independent 
of the size of the simulated system. The computational scheme allows for efficient calculations of 
the superfluid fraction and off-diagonal correlations in space-time, for system sizes which are orders 
of magnitude larger than those accessible to conventional PIMC. We present illustrative results for 
the superfluid transition in bulk liquid ''He in two and three dimensions, as well as the calculation 
of the chemical potential of hep *He. 
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I. INTRODUCTION 

It is now twenty years since Ceperley and Pollock (CP) 
carried out the first Path Integral Monte Carlo (PIMC) 
simulation of the superfluid transition of liquid ^Hei Al- 
beit restricted to a system of 64 '^He atoms with periodic 
boundary conditions, that study demonstrated the fea- 
sibility of ah initio numerical studies of quantum many- 
body systems, the mass of the particles and the interac- 
tion potential being the sole input to the calculation. 

The PIMC method, in the form developed by CP 
(henceforth referred to as "conventional"), has since 
played a major role in the theoretical investigation of 
quantum many-body systems. Not only has it provided 
quantitative results for a wide range of physical systems, 
it has also shaped, to some extent, our qualitative un- 
derstanding of such phenomena as superfluidity (SF) and 
Bose condensation, at the microscopic level. At least for 
Bose systems, PIMC is the only presently known method 
capable of furnishing in principle exact numerical esti- 
mates of physical observables at finite temperature (T) , 
including the superfluid (ps) and condensate (uo) frac- 
tions. Moreover, despite the notorious sign problem, 
that has so far made it impossible to obtain equally high 
quality results for Fermi systems, PIMC proves a valid 
option in this case as well, allowing one to obtain approx- 
imate estimates of accuracy at least comparable to that 
afforded by the other leading methodsiSi^ 

It thus seems reasonable to regard PIMC as a realistic 
option to investigate ever more complex quantum many- 
body systems, and it makes sense to try and overcome 
its most important present limitations. Aside from the 
above-mentioned sign problem, which we do not discuss 



in this paper, the main bottleneck of the current PIMC 
technology is inarguably the maximum system size (i.e., 
number N of particles) for which accurate estimates can 
be obtained, in a reasonable amount of computer time. 
Specifically, the computational effort required to study 
properties that most directly depend on particle indis- 
tinguishability, is observed to scale prohibitively with N . 

For example, the superfluid fraction ps is obtained in 
a PIMC simulation of bulk condensed matter, by means 
of the so-called winding number estimator,^ which can 
only take on a nonzero value if long permutation cycles 
of identical particles occur. In conventional PIMC, the 
frequency with which such cycles are sampled, is an ex- 
ponentially decreasing function of N . For this reason, 
and in spite of (at least) a hundredfold increase in com- 
puter speed)^ since the pioneering work of Ref. Q it has 
not proven possible to obtain estimates of the superfluid 
fraction in bulk liquid ^He for finite systems of more than 
N—QA particles. Besides the unfavorable scaling of com- 
puting resources as a function of N , another major is- 
sue that this entails is the difficulty of assessing reliably 
whether the observed absence of long permutation cy- 
cles refiects a genuine physical effect, or merely lack of 
ergodicity of the path sampling scheme.^ 

How important is the above size limitation ? For 
most observables diagonal in the coordinate represen- 
tation, one can often approach surprisingly closely the 
thermodynamic limit by simulating as few as ~ 30 parti- 
cles, especially for systems characterized by short-ranged, 
Lennard-Jones-type interactions. On the other hand, 
an accurate quantitative characterization of the super- 
fiuid transition (including the calculation of the transi- 
tion temperature T^) can only be obtained via finite-size 
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scaling analysis of results for Ps{T) and/or no(T). The 
reliability of this procedure crucially hinges on the avail- 
ability of data for large systems of significantly different 
sizes. Attempts to estimate for supcrfluid ^He, based 
on PIMC data for Ps{T) for systems of size iV=64 and 
smaller, failed to yield quantitative results^ 

But there are other reasons, arguably more important 
than the mere pursuit of numerical accuracy, pointing to 
the importance and timeliness of extending by one or two 
orders of magnitude the size of the systems accessible to 
PIMC. Quite generally, in order for a scientific question 
to be meaningfully addressed by numerical simulations, 
the size of the simulated system should be greater than 
all characteristic length scales affecting the physics of in- 
terest. This is particularly important in the study of in- 
homogeneous phases of matter, or fluids in confinement, 
or restricted geometries. An example is provided by the 
study of helium fluid in porous glass, such as Vycor; the 
diameter of a characteristic pore is of the order of a few 
tens of A. Thus, a realistic calculation, at the typical liq- 
uid helium density, requires that one be able to simulate 
a system comprising several thousands of atoms. Other 
examples are the numerical investigations of multicompo- 
nent systems, as well as of grain boundaries, dislocations 
and other defects in quantum solids, or of incommensu- 
rate phases of films of helium or para-hydrogen adsorbed 
over substrates such as graphite. 

Over the past two decades, there has been relatively 
little experimentation with approaches to PIMC simula- 
tions differing in some important aspects from the con- 
ventional one of CP, thoroughly described in Ref. 0- 
As mentioned above, in conventional PIMC the simu- 
lation of properties that are most directly affected by 
quantum statistics (i.e., by many-particle permutations), 
suffers from a very unfavorable scaling of required com- 
puter time with system size. This hurdle seems difficult 
to conquer within conventional PIMC, and more gener- 
ally within any Monte Carlo scheme formulated in the 
canonical ensemble, in which the winding number be- 
comes "topologically locked" in the N ^ oo limit. ^ 

On the other hand, the same hurdle has been com- 
pletely overcome in Quantum Monte Carlo (QMC) simu- 
lations of lattice models. A lattice Path Integral scheme 
based on an alternative sampling approach, known as 
worm algorithm ( WA) ^ has been demonstrated to allow 
for efficient calculations of winding numbers and of the 
one-particle Green function G, for systems of as many as 
^ 10^ particles.^ It is particularly useful for the studies of 
critical phenomena since it does not suffer from the criti- 
cal slowing down problemifi present in other local- update 
schemes. 

The WA has been recently extended to the study of 
systems in continuous space; it has first been shown 
to afford the simulation of the superfluid transition in 
liquid ^He in two dimensions, for systems comprising as 
many as 2500 particles, i.e., about 100 times greater than 
those accessible to conventional PIMC. Subsequently, it 
has been applied to the study of Bose condensation in 



crystalline ''He,— as well as to the investigation of super- 
fluid properties of para-Hydrogen droplets 

In all of these applications, the WA has provided ac- 
curate numerical results, simply not obtainable with any 
other existing method. It need be stressed, however, that 
the WA is not merely about doing large system sizes 
(important as this is); it is also the first grand canon- 
ical QMC method with local updates to incorporate in 
full quantum statistics. It affords the exact computation 
of imaginary-time off-diagonal correlations, such as the 
onc-particlc Matsubara Green function, that are not ac- 
cessible to conventional PIMC (nor to any other QMC 
technique in continuous space). 

In this paper, we provide a detailed description of 
this new, powerful computational tool, which promises 
to open novel avenues to the theoretical exploration of 
strongly correlated many-body system. The manuscript 
is organized as follows: in the next section iQlJ, we de- 
scribe the simplest implementation of the WA, also intro- 
ducing our nomenclature, configurational space and data 
structure. In section Hill we discuss in detail the compu- 
tation of the Matsubara Green function, as well as of its 
equal-time limit, namely the one-particle density matrix. 
In section IIVI we describe an important enhancement 
of the simple implementation, namely the expansion of 
(the attractive part of) the potential energy of interaction 
among particles in elementary (diagrammatic) contribu- 
tions, that are then sampled by a Monte Carlo method. 
This scheme, which falls in the general category of Dia- 
grammatic Monte Carlo techniques,^"* allows for the full 
inclusion of the contribution of the potential energy (for 
an important class of potentials) at a computational cost 
that is independent of the size of the system. 

In Sec. we provide some quantitative information 
related to the specific utilization of the WA in simulation 
studies of superfluid ^He. In Sec. I VII we offer a quantita- 
tive demonstration of the power of the WA, by illustrat- 
ing in detail our results for simulations of the superfluid 
transition of liquid '''He in two and three dimensions. We 
outline our conclusions, and discuss outlook for future 
applications of the WA, in Sec. IVIII 



II. SIMPLEST VERSION 

We begin with some basic notation. We assume for 
deflniteness a system of identical particles (in d dimen- 
sions) obeying Bose statisticsii^ Let m be the mass of 
each particle. The system is enclosed in a cubic vessel 
of volume V = L'^, with periodic boundary conditions 
in all directions (other geometries, boundary conditions, 
and / or external forces require obvious and minimal mod- 
iflcations which are standard for any QMC scheme). 

We make from the outset the assumption of working in 
the grand canonical ensemble; that is, the system is held 
in thermal equilibrium with a heat reservoir at tempera- 
ture T — 1/ (3 (we set fcB = l), with which it can exchange 
particles as well. Consequently, in order to specify the 



3 



thermodynamic state of the system, we need to assign 
the chemical potential ^, which is an input parameter in 
our computational scheme, just like the temperature T. 
The number of particles N is allowed to fluctuate. 

Let H be the (many-body) system Hamiltonian, which 
we assume of the following form: 



N 



(2.1) 



where A = ?i^/2m and where is a pairwise interaction 
potential that depends only of the relative distance be- 
tween any two particles>i& Below, R will always be used 
as a collective "coordinate" , representing positions of all 
particles in the system, i.e., R = (ri,r2, . . . 



A. Configurational Space 

A fundamental aspect of the WA, which crucially dis- 
tinguishes it from conventional PIMC and from all ex- 
isting QMC methods in the continuum, is that it op- 
erates in an extended configurational space, containing 
both closed world line configurations (henceforth referred 
to as Z- or diagonal configurations), as well as config- 
urations containing one open world line (worm). The 
Z-configurations contribute to the partition function, 
whereas those with an open world line contribute to the 
one-particle Matsubara Green function; in the following, 
the latter will be referred to as G- (or, off-diagonal) con- 
figurations. As we shall see, all topologically non-trivial 
modifications of world lines occur in the off-diagonal con- 
figurational space (or, G-sector). The sampling process 
allows for transitions from the G- to the Z-sector (by 
closing, or removing the existing open world line) and 
vice versa (by creating a new open world line, or by 
opening an existing closed one). Expectation values of 
all physical quantities of interest (with the exception of 
the Green function), including particle and winding num- 
bers, are only updated when the random walk generates 
a diagonal configuration. 

Next, we proceed to describe in detail the two sectors 
in which our configuration space is conceptually divided. 



1. The Z -sector 

The Z-sector of our configuration space, is nothing 
but the full configuration space of conventional PIMC. 
It naturally emerges from the Path Integral representa- 
tion of the grand partition function Z = Tre"'"^, where 
k = H - ^iN. 

Each Z — configuration is a discrete imaginary-time 
many-particle path, X = {Rq, Ri, R2, . . . , Rp), with 
Rp = Rq (except for a possible permutation of particle 
labels), representing the integrand in the asymptotically 



exact (in the P — > 00 limit) integral decomposition of Z 
Z^Y. e^""^ / '^^ ^(^' ^) e"''^"'^ ' (2.2) 

iV=0 

where dX = dRodRi...dRp^i, e — P/P, and where 
p-i 

A{X,e)= l[pF{Rj,R,+,,e) . 

j=o 

In turn, pp is a product of free-particle imaginary-time 
propagators, i.e., with an obvious notation, 

N 

PF{Rj,Rj+i,e) = ]J/9o(ry,ri,j+i,e) (2.3) 



with 



Po(r,r',e) = (47rAe)-'^/2 



exp 



(r-r') 



l\2 
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(2.4) 



The function U in Eq. H2.2|l incorporates correlations, 
both in space and in imaginary time, arising from inter- 
actions among particles. U is chosen so that, in the £ — > 
limit, the distribution of discrete paths X will asymp- 
totically approach the correct continuous limit. Several 
choices are possibleSi for U, but the simplest version of 
the algorithm described in this section does not depend 
on its particular form. In what follows, we refer to the 
product W{X) = A{X,e) e"'^^'^^ as a configurational 
weight. 

Eq. H2.2|l implies the following configuration space 
structure: one has N single-particle paths (world lines), 
labeled i — 1,2, ...,7V, propagating in the discretized 
imaginary time interval [0,/3] (specifically, to = 0, tp = 
(3). Each world line consists of P successively linked 
"beads" (particle positions), labeled by the index of the 
corresponding imaginary time "slice" , j = 0,...,P — 1. 
The j-th bead of the z-th world line is positioned at . 

As a result of /3-periodicity, coupled with the physi- 
cal indistinguishability of particles, the (P — l)-st bead 
of each world line must be linked to the zero-th bead of 
either the same, or another world line. For both theoret- 
ical and practical (data structure) purposes, it is advan- 
tageous to guarantee /3-periodicity automatically; to this 
aim, we regard world line configurations as closed loops 
on a ((i+ l)-dimensional surface of a (d-l- 2)-dimensional 
/3-cylinder, on which lie P equidistant (and equivalent) 
imaginary time "hyperplanes" (corresponding to the dif- 
ferent time slices), labelled j = 0, . . . , P — 1. It should be 
noted that the total number of world line loops defined 
on the /3-cylinder can be different from the total num- 
ber of particles; this is because a single world line which 
"winds around" the imaginary time interval / times be- 
fore returning to its initial position represents not just 
one particle, but rather I particles involved in the same 
exchange cycle. The presence of such exchange cycles is 
essential, in order to incorporate in the computational 
scheme the symmetry of the system with respect to par- 
ticle permutations. 
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2. The G-sector 

The G-sector of our configurational space comes from 
the representation — analogous to that of the partition 
function, Eq. (|2.2|) — of the one-particle Matsubara Green 
function 

.g(ri,r2,T) 



Z 



G(ri,r2,r) = {T{^{v,,t) i^Hr^^O)}) = 

' (2.5) 

where (. . .) denotes thermal averaging, T is the time- 
ordering operator and ■\\^(x:^t) and ■!/'(r,r) are (Bose) 
particle creation and annihilation operators in Matsubara 
representation. The structure of the integral representa- 
tion oig(y\ , r2 , r) is very similar to that of Z. In fact, the 
only qualitative difference of a G-sector ("off-diagonal") 
configuration from a diagonal one, is that the former con- 
tains a worm, that is, a world line on a /3-cylinder with 
two ends — the "head" and the "tail" — corresponding to 
the Green function annihilation and creation operators, 
respectively. The two special beads at the open world line 
ends are named (for historical reasons) Ira (T) and Masha 
(Ai). Configurations in which X and A4 are located in 
space-time at points (ri,Ti) and (ymtTm) contribute to 
g{rx, r_A4, tj — tm) with the weight defined in accordance 
with Eq. (|2.2I) generalized to include the off-diagonal con- 
figuration sector. 

Formally, the ensemble of WA configurations corre- 
sponds to the generalized partition function 



Z' 



(2.6) 



with 



Z' = dn dvM g{rx, Tm, e(ji - Jm)) ■ (2.7) 



3x,3M 

The value of dimensionless parameter G only affects the 
efficiency of the simulation, as G controls the relative 
statistics of Z- and G-sectors; for the moment, we leave 
it undetermined, to come back to it later on, when dis- 
cussing updates. 

An important new feature that arises when going from 
Z to Z\Y, is that the number of continuous variables in 
G-sector configurations is not constant, but rather varies 
from configuration to configuration. This immediately 
points to Diagrammatic Monte Carlo^^ as a general way 
to perform updates whenever the number of variables to 
sample, is itself variable. 

The sampling of paths {X;}, is implemented within 
the WA exclusively through a set of simple, local space- 
time updates involving X (or, A4). The particle number 
becomes configuration- and time-dependent (there is one 
less particle between X and A4, than in the rest of the 
path). This clearly shows how, by its very construction, 
the WA opens up the possibility of working in the grand 
canonical ensemble, with the chemical potential fj, being 
an input parameter»ii Obviously, WA updates can be 
combined with the conventional PIMC updates to have 
the most flexible scheme. 



In accordance with Eq. (|2.7() , the simplest estimators in 
WA are [assuming that configurations are sampled from 
the probability density VF(-'^)] 



^(2) _ I 1 , if in Z-sector , 
, if in G-sector , 



j(G) — J ' Z-sector , 

1 , if in G-sector . 



(2.8) 



(2.9) 



In the statistical limit, their Monte Carlo averages are 
('5'''')mc = Z/Zw , (2.10) 



MC 



CP 



■jW 



^ Mridr2 5(ri,r2,ej) , (2.11) 

A n J 



In particular 



((5(^)) 



MC 



CP Y /rfridr2 5(ri,r2,£j) . (2.12) 



The simplest estimator for g(r,T), r = ej, is given by 

i^^^^ ^j,Ui-JM)'^i'^'i- -i"i)(5(r2 -rA4))MC = 
CP 

= _g(ri,r2,ej) . (2.13) 

However, below we introduce a more elaborate scheme 
which allows one to circumvent the problem of working 
with generalized functions (typically solved by collecting 
statistics to finite-size spatial bins at the expense of an 
additional systematic error). 



B. Data Structure and Updates 

In this section, we describe a set of ergodic local up- 
dates which sample the extended configuration space, 
switching between the Z- and G-sectors. Updates which 
change the number of continuous variables in X, are ar- 
ranged in complementary pairs, designed so as to satisfy 
the requirement of detailed balance. General principles of 
balancing complementary pairs can be found in Ref. IT^ 
We have three pairs of updates altogether: Open/Close, 
Insert/Remove, and Advance/Recede. Only the Swap up- 
date in the list below does not fall in this category, be- 
cause it preserves the number of variables, i.e., it is self- 
complementary. Proposed updates are either accepted or 
rejected based on the Metropolis algorithm, according to 
the standard procedure,^^ 

In order to keep the presentation simple, the sampling 
scheme described below is one in which every update can 
be proposed, regardless of its apphcability to the cur- 
rent configuration (for example, the proposal to Remove 
the worm is allowed even if there is no worm in the cur- 
rent configuration, in which case the proposed update will 
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FIG. 1: World line beads on the /3-cylinder and nearest neigh- 
bor associations between them. 



necessarily be rejected). It should be understood, how- 
ever, that standard sampling tricks can be used, whereby 
only applicable updates are proposed, thereby enhancing 
the performance«i^ 

To be specific, let us adopt the following data struc- 
ture: all beads are labeled, and each bead is linked to 
its two world-line neighbors, the next and the previous 
beads on the /3-cylinder, see Fig. ^ It proves convenient 
to introduce two functions, 'nexf and ^prev\ mapping 
each bead onto its next and previous neighbor, respec- 
tively. Correspondingly, a — next"^{a) means that the 
bead a is the result of the m-fold application of 'nexf to 
the bead a; likewise, a = prev"^{a). 

In order to have an efficient spatial addressing of beads, 
and thus a scalable algorithm (i.e., one in which the num- 
ber of operations required to perform updates does not 
depend on either system size, nor temperature), we use 
nearest- neighbor tables at all imaginary time slices ,2fl, 
That is, for each time slice, the volume of the system 
is hashed into equal microscopic "bins", labeled by the 
discrete time label, j = 0, 1, 2, . . . , (P — 1), and discrete 
radius vector, TZ. [In practice, we set the bin volume, Q, 
to be of the order of the volume per particle.] For each 
bin there is a list of all beads contained in it. Accord- 
ingly, addressing/searching relevant beads is performed 
through nearest-neighbor tables and next /previous links. 

(la) Open. The update is only possible if the configu- 
ration is diagonal. A bead a is selected at random. An in- 
teger number M is selected at random within the interval 
[1, M] with M < P being an arbitrary algorithm param- 
eter. Then, {M — 1) beads, namely, next^{a), next'^{a), 
. . . , nexi^*^^^-' (a) are removed, so that a worm appears 



with J — a and ^ a = next 
probability for this update is 



(a). The acceptance 



AU-iiMe 



P =min|l ^^^bd e 



(2.14) 



where At/ = U{X) — U{X*) is the difference between 
the {/-function values for the initial (X) and proposed 



{X*) configurations, iVbd is the total number of beads in 
the initial diagonal configuration equal to the number of 
particles, N, times the number of slices, iVbd = NP. If 
M — 1, then no beads are removed; only the link between 
the beads a and a disappears with the appearance of I 
and M. 

(lb) Close. This update is only possible if the con- 
figuration is off-diagonal. Let the integer M > be 
the discrete algebraic distance from X to A^, understood 
as a number of time-slice steps — in the positive direc- 
tion on the /3-cylinder — one has to make to reach the 
slice at which Ai is currently positioned, starting from 
that of J. If M > M or M = 0, then the move is 
rejectedjiS otherwise, one proposes to generate a piece 
of world line connecting X to A^, thereby rendering the 
configuration diagonal. If M > 1, the corresponding spa- 
tial positions of new (M — 1) beads, ri, r2 . . . , tm-i, are 
sampled from the product of M free-particle propagators 
ni!liPo(ri._i,r^,e), where Tq = rj and ym = tm- The 
probability of accepting the move is 



Pel = min<^ 1, 



Po(rj,r^,Me) e' 



AU+fj.Me 



CMNhd 



(2.15) 



where Nm is the number of beads in the final diagonal 
configuration. In our implementation, proposed Open 
and Close updates are automatically rejected whenever 
the quantity 

(ri - tm^ 



4AfAe 

is larger than some (arbitrary) number of order unity, so 
as to avoid small acceptance ratios in the close update 
when the worm ends are far away in space (in our simu- 
lations we set this number equal to 4). 

(2a) Insert. The other way to create an off-diagonal 
configuration from a diagonal one, besides Open, is to 
seed a new, M-link long, open world line. The number 
of links 1 < M < M and the position of A4 in space- 
time are selected at random. The spatial positions of 
the other M beads are generated from the product of 
M free-particle propagators. The move is accepted with 
probability 



Pi„= min{l, CYPMe^^+f'^^''} 



(2.16) 



where V is the volume of the system, as mentioned above. 

(2b) Remove. The removal of the worm, i.e., of the 
world line connecting A4 to I, is proposed, provided its 
algebraic length is 1 < AI < M. (If M > M, the proposal 
is rejected^) The acceptance probability for the move is 



Pim = min{l, e 



/CVPM} 



(2.17) 



At this point, we are in position to discuss the value of the 
constant C, which up to now we have left undetermined. 
A natural choice is 



C^Cq/VPM, Co ^ 0(1), 



(2.18) 



so that the probabiHties to Open, Close, Insert, or Re- 
move a worm, do not contain macroscopically large/small 
factors, and are of order unity at the optimal choice of 
M. Normally, optimal M is such that the time eM is of 
the order of the characteristic single-particle time, which 
guarantees that the exponentials in the acceptance prob- 
abilities are of order unity, while the propagators are of 
the order of the particle number density. This implies 
the following scaling: 



(3a) Advance. This move advances X by a random 
number M of slices forward in time. Its implementation 
is similar to that of Insert. The acceptance probability 
is 

Pad = min{l, c^u+^.Me^^ (2.20) 

(3b) Recede. Now X is displaced backwards in time 
(again, in a /3-periodic sense), by erasing M consecutive 
links; the number 1 < M < M is selected at random. 
The acceptance probability is 

P,c = min{l, e-^^-^*^^} (2.21) 

If M turns out to be equal to or larger than the number 
of links in the worm, the update is rejectediiS 

(4) Swap. This update is applicable to off-diagonal 
configurations only and is illustrated in Fig. |21 Let T be 
positioned on the j-th slice in the bin (7?i, j). Consider 
the {j + M)-th slice (because of /3-periodicity, this ad- 
dition is understood modulo P) and create a temporary 
list, Cx, of all the beads that are contained, at the slice 
j + M, in the bins that spatially coincide with the bin 
{T^XtJ) or with one of its nearest neighbors. Select one 
bead, a, from the list, with the probability 

Ta = Po(ri,r„,Me)/I]i , (2.22) 

where 

Si - ^ Po{Yx,Y,,Me) (2.23) 

is the normalization factor. If any of the beads a, 
prev^{a), prev'^{a), . . ., prev^\a) coincides with M, the 
move is rejected. Next, consider the bead C, = prev^^ {a) 
and identify its bin, {TZQ,j). One must now check 
whether the bead a is contained in the bin spatially co- 
inciding with (TZQ,j), or with one of its nearest neigh- 
bors; if that is not the case, the move is rejected. A 
second list is then created, Cq, of beads contained in 

the {jtc,ij + ^) bin and in its nearest neighboring ones. 
At this point, a set of new positions between beads a 
and I is generated (in the same way as in the Close 
move), Ti, . . . , rjg_]^ for the beads prev^{a), prev'^{a), 
. . ., prev^*^""'^-' (a), respectively. One may now rename I 
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FIG. 2: An illustration of the Swap update and the notation 
introduced in the main text. 

into and vice versa, and re-link beads in the following 
manner: former prew*-*^"'"^-' (a) becomes prew (I), etc. As 
a result of this world-line reconnection, a piece of world 
line between the original bead C and the bead a is erased, 
while a new piece of world line appears that connects the 
former bead X with the bead a. The move is accepted 
with probability 

Psw - min{l, e'^^Ei/Ej . (2.24) 

Here T,q is defined similarly to Sj in Eq. (|2.23() using the 
£^-list. 

The Swap move generates all possible many-body per- 
mutations through a chain of local single-particle up- 
dates. Since no two particles need be brought within 
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a distance of the order of the hard core of a typical in- 
teratomic potential potential, this move enjoys a high ac- 
ceptance rate, similar to that for the Advance/Recede up- 
dates (we provide a quantitative example of acceptance 
ratios measured in simulations of the two-dimensional 
^He liquid below). It must be emphasized that in our 
algorithm, unlike in conventional PIMC, arbitrary per- 
mutations of identical particles, as well as macroscopic 
exchange cycles need not be explicitly sampled. For, they 
appear automatically, if the physical conditions warrant 
them. This is because the statistics of the relative po- 
sitions for the worm ends is given exactly by the Green 
function G(ri,r2,T). 

In a typical sweep, a worm is created (either by opening 
an existing closed world line, or by inserting a new one), 
advances and/or recedes in imaginary time and performs 
a number of swaps, until it finally closes or is removed. It 
is easy to convince oneself that, as a result of this proce- 
dure, an exchange cycle involving a macroscopic number 
of particles can appear in just one sweep. 

Since the complementary pairs are detail-balanced, the 
final results do not depend on the global probability of 
addressing each pair, so long as the probabilities of ad- 
dressing each update within a complementary pair are 
equal, which is assumed in the acceptance probabilities 
presented above. Otherwise, the acceptance probabilities 
should be modified as follows: if the probability to ad- 
dress update A is ua and the probability to address a 
complementary update is ub then Pa Pa{ub /ua)- 

As noted in the Introduction, the description of WA 
given above is complete even without expansion of the 
pairwise potential tail into diagrams (described below). 
In the original version of the code we did not use this effi- 
ciency enhancing modification of the conventional PIMC 
configuration space, and still were able to perform accu- 
rate studies of several hundred atoms. Since expansion 
into diagrams works only for the attractive part of the 
potential it can not be used for purely repulsive models. 



III. PHYSICAL ESTIMATORS 

The statistical estimators for all physical quantities 
that are computed in the Z-sector, including those of 
all energetic and structural properties, are identical with 
those utilized in conventional PIMC. We therefore refer 
the reader to Ref . for a detailed discussion of these esti- 
mators. Instead, we focus our discussion here on the esti- 
mator that is a characteristic feature of the WA, namely 
that of the Matsubara Green function. 



A. Green function and density matrix estimators 

Since the Green function is sensitive only to relative 
distances in time, without loss of generality we can fix 
JM = to simplify the notation. For an extra simplic- 
ity, we will also assume that at least one of the following 



two statements is true: (i) The problem is translation- 
ally invariant in the coordinate space, (ii) The quantity 
of interest is the Green function averaged over spatial 
translations, 

G{v,t) = j dr'G(r + r',r',r) . (3.1) 

In both cases we can formally fix tm = 0. [The general- 
ization of the treatment to the case of two independent 
spatial coordinates is straightforward.] 

Consider the configurations, or path-integral diagrams, 
'D^irx, ji) contributing to the function g{r,ej): 

g{vi,n) = ^Pe(ri,ji). (3.2) 

e 

The subscript ^ stands for all variables of the diagram 
except for the end point space-time positions. That is ^ 
contains both positions of the beads (in which case the 
summation is understood as integration) and a topolog- 
ical structure of the world lines on the /3-cylinder. 

Comparing the diagrams with different jx and rj (in- 
cidentally, it is precisely this comparison that stands be- 
hind the acceptance probabilities of Advance update), we 
readily see that given some diagram T>(^^{tq, j^) we can 
"upgrade" it to a diagram 'D^{v,j) with j > jo by at- 
taching to Z a world line piece of M = j — jo beads. 
Specifically, 

T^d^J) = ^'Coli'o, jo) i?joj(ro,r) (jo < j) , (3.3) 

M 

i?,„,(ro,r) = e^^+^-*^ [] Po(r.^i,r„e) , (3.4) 
i/=i 

where ^ = {^o, ri, r2, . . . , tm-i}, tm = i", and the mean- 
ing of AU is the same as in Advance update: AU — 
U^o ^ U^, where and correspond to the diagrams 
and -D^qi respectively. Hence, we have 

.9(1", ej) = X! dro... dvM-i V^„{ro,jo) Rj„j{ro,r) = 

Co 

Jfj^j d^o ■ ■ - drM-iV^oiT^oJo) Rjoji^o,r)S^j^j\ (3.5) 

5o Jo 

where 

c(M) ^ / 1 , if jo e [j - M, j - 1] , , . 

\ , otherwise . ^ ' 

Though Eq. (|3.6(l is not the final answer yet, it al- 
ready contains an important element. It allows one to 
sample the time slice j from the adjacent time slices 
jo G [j - M,j -1]. The problem with Eq. (|^ is that 
we still have a continuous variable r while we would like 
to know the Green function only at the discrete set of 
pre-defined points {fp}. To proceed further, we utilize 
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(and to a certain extent generalize) the idea which has 
been already used in diagrammatic Monte Carlo»^ Sup- 
pose we are interested in g(rp,ej), where Tp G is a 
pre-selected point for collecting statistics and is a 3D 
volume containing this point. We formally rewrite (|3.5() 
at point Tp as (below tm = fp) 



io do 



(M) 

ioi 



X / drWMiro,r)S; 



(Vp) 



X / rfri . . . dru-i , ttt Q 

Po(ro,rp,£M) 



(3.7) 



Q 



(MV Po(''Oi''pig^^) AtZ+MeA/ o\ 




if r e , 
otherwise . 



(3.9) 



The value of AC/ corresponds to the extra piece of world 
line (ro, ri, . . . , tm-i, rp). In principle, Wm{yq, r) is an 
arbitrary function, but we want it to be positive-definite 
and normalized (the integration is over the whole system 
volume, not just Vp), i.e. to have the meaning of the 
probability density 



j dvWM{Y0,Y) = 1 



(3.10) 



Moreover, for our purposes the most reasonable choice is 
simply 



Wj\/(ro,r) = Po(ro,r,eM). 



(3.11) 



Now we just need to interpret the relation H3.7(l in terms 
of a stochastic process. First, we integrate/sum this re- 
lation over the position of Ai in space/time and com- 
pensate for that by dividing H3.7(l by PV . To interpret 
the first line in (|3.7|) . we recall that in accordance with 
(|2.7() the value of is equal to the probability den- 
sity to sample the corresponding diagram of the G-sector 
times the factor Zyt/jC. Hence, we can interpret the 
first line as averaging — over the ensemble of all Monte 
Carlo diagrams — of the stochastic variable given by the 
rest of the expression times projector S^^^ times projector 

"^fe-jAi) j ^i™'^^ I C^P- The second line says that the 
evaluation of the stochastic variable starts with sampling 
a vector r distributed in accordance with Eq. H3.11|l . The 
projector b^^^ means that if r ^ "(^, then the stochastic 
variable is automatically zero. If r G T^, then the eval- 
uation procedure continues in accordance with the third 
line of Eq. (|3.7|) . which we interpret as sampling (Af — 1) 
auxiliary variables, ri, . . . , Tm-i, where M = j —ji+jM- 
(In the special case of M = 1 , auxiliary variables are not 



FIG. 3: Construction of the estimator for the Green function 
G(rp,ej), where without loss of generality we assume that 
rM = 0, Jm = 0. 



sampled.) The auxiliary variables are sampled from the 
distribution 



Po(ro,ri,£)po(ri,r2,£) . . .po(rM-i,rp,g) 
Po(ro,rp,eM) 



(3.12) 



When these are fixed, the stochastic variable in 
question — up to the global pre-factors discussed above — 
is nothing other than Q, defined by Eq. I|3.8(l . 
Summarizing, we have derived the estimator 



CVP 



(3.13) 

where the variable Q is calculated in accordance with 
Eq. H3.8() in terms of the auxiliary variables of the above- 
described procedure, which we illustrate in Fig. |31 

Q ^ } po{ra,rp,eM) au+^,,m 
MVp po{ro,r,eM) 

The free parameters of the procedure, M and Vp, are 
optimized to yield the best possible convergence. 

Our final note is that if one is interested in the Green 
function at a certain momentum p then the correspond- 
ing estimator does not require any elaboration presented 
above due to extra integration over r. Indeed, according 
to Eq. (|2.13l) we have 



MC 



CVP 
Zw 



dip, £j) • 



(3.14) 
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FIG. 4: Configuration space witii diagrammatic bonds be- 
tween the equal-time beads. 



IV. ENHANCED VERSION: DIAGRAMMATIC 
EXPANSION OF THE ATTRACTIVE 
POTENTIAL TAIL 

The WA described in^Jcan be (and, has been) used for 
efficient simulations of systems comprising a few hundred 
particles. This section describes a general procedure that 
significantly improves performance of PIMC, through an 
efficient sampling of the contribution of the attractive 
tail of the pair potential. The trick per se is not directly 
related to the WA, and can be implemented within other 
PIMC schemes. Obviously, the modification of the con- 
figurational space it brings about should be adequately 
taken into account in the updates. Moreover, in order for 
the trick to be applicable an attractive tail of the inter- 
particle potential must exist, i.e., what described below 
does not apply to purely repulsive interactions, such as 
hard-sphere. 

Straightforward schemes for calculating the potential 
energy exponent of a multi-particle path, U{X), have 
to deal with a compromise between accuracy and per- 
formance. By a "straightforward scheme" , we mean one 
that treats the pair interaction between any two particles 
on the same footing, irrespectively of the value of the po- 
tential, as long as it is non-zero, which is the case in all 
realistic models. The computational complexity of such 
schemes scales linearly with N per single-particle update. 
To improve on efficiency, one typically truncates the po- 
tential at some distance, thus introducing a systematic 
error. 

There is, however, a way to reduce radically the com- 
putational effort required for accurate treatment of the 
potential tail, without introducing any additional system- 
atic error. The only price to pay is a more complex con- 
figuration space (see Fig. ^ and, correspondingly, addi- 
tional updates to sample it. 



We confine ourselves to the simplest case, i.e., the one 
in which U{X) is a sum of pairwise contributions: 



p-i 

3=0 {a,b,) 



(4.1) 



where beads on the j-th time slice are labeled using sub- 
scripts aj and bj, the notation (ajbj) stands for all pairs 
of beads on a given slice, and r's are the spatial co- 
ordinates of the beads. In the simplest^ choice for U 
u{r) = v{r)e. Correspondingly, 



e-^ = 



nn e 

3 {ajhj) 



■u{r 



(4.2) 



An important observation now is that if |ra. — rf,.| > 
Tc, where is some distance, greater than the radius 
of the repulsive core of the pair potential, then uiva- — 
Vb- ) < 0, and the corresponding pairwise exponential can 
be identically split into a sum of two positive definite 
terms: 



-ti(r„ -Ti, 



1 



1 



(4.3) 



Graphically, this decomposition can be represented as fol- 
lows: Any two beads within one and the same slice with 
I Taj — ffcj I > f'c, now may or not share an extra graphical 
element, a 6onc?, see Fig.0] The new configuration space 
is reminiscent of the Feynman's diagrammatic expansion 
in powers of the interaction potential. The absence of a 
bond between two beads Oj and hj represents the first 
term in Eq. i.e. unity, and means that the beads 

at a distance |ra. > Tc do not interact. A bond 

between beads aj and bj represents the second term of 
Eq. 1)4.3(1 . Its relative contribution to the statistics as 
compared to the case of no bond has an extra factor of 



(bond factor) 



(4.4) 



If one of the beads in the bond is a worm, then the bond 
factor is 



-ii(r„ -ri,.)/2 1 
e ^ 1 "j" — \ 



(worm bond factor) . (4.5) 



Formally, we attribute potential energy exponents to 
beads, but in reality they represent interactions between 
the world line trajectories; since the worm bead has a 
trajectory attached to it from one side only, its potential 
energy exponent is reduced by a factor of two. Corre- 
spondingly, the two worm beads do not interact and thus 
T can not be connected by the bond to Ai. 

The standard MC prescription for sampling the new 
configuration space would be to have updates which cre- 
ate and remove bonds. It is easy to see how a "rad- 
ical gain in performance" can be achieved, as, statis- 
tically, the probability for two beads aj and bj with 
\ra — ffo I > Tc to share a bond is much smaller than unity. 
Indeed, this probability is proportional to the bond factor 
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(|4.4() . which can be estimated as ~ v{r)e <C 1. In order 
to sample the new configuration space, we introduce a 
pair of complementary updates that create and remove 
bonds. This pair of updates is reminiscent of diagram- 
matic Monte Carlo updates that create and delete beads. 
There is, however, an important mathematical difference. 
When creating new beads one has to seed new continu- 
ous variables associated with them. Bonds are created 
between existing beads, and formally the new updates 
are of the standard Metropolis type. 

Create Bond. The update is only possible if the config- 
uration is off-diagonal. An integer number M is selected 
at random within the interval [0, M]; M < P. This num- 
ber is used to select the first bead, aj = prev^-' {I), in 
the pair to be connected by the bond. If prev^^ (J) is not 
defined because the world line terminates at M, the up- 
date is rejected. The second bead, bj, is selected in two 
steps. First, within the slice j of the bead Uj we select a 
spatial bin, B. This is done by randomly generating the 
bin label from some probability distribution, P^g, which, 
in general, depends on the distance between B and the 
bin A that contains the first bead aj . (We discuss a rea- 
sonable choice of P_ab below.) Let be the number of 
beads in the bin B. If rig = 0, the update is rejected. 
If ng > 0, we select at random a bead from the bin B 
and call it bj . If it happens that aj = X and bj = M. the 
update is rejected because physically the world line ends 
do not interact. Also, if Oj and bj are already connected 
by a bond, or the distance between the selected beads is 
smaller than Tc, the update is rejected. Otherwise, we 
propose to create a bond between the selected pair of 
beads, and accept the proposal with the probability 



-Pcrb - 



{M + \)nB 



bnd 



1)Pab 



-/"(raj-rtj) _ ^ 



(4.6) 



Here ^bnd is the total number of bonds in the initial 
configuration associated with the beads I, prev{I), . . ., 
prev^\l) or X, previl), . . ., if prev^^T) is not de- 
fined. The balancing factor (/bnd + 1) naturally emerges 
from the Remove Bond update which is complementary 
to the Create Bond. An additional factor / = 1/2 in the 
exponent is necessary only if one of the beads in the pair 
is the world line end; otherwise, f — 1. 

Remove Bond. The update is only possible if the con- 
figuration is off-diagonal. We list all l^nd bonds asso- 
ciated with the beads X, prev{X), . . ., prev^ (X) or X, 
prev{X), . . ., M ii prev^{X) is not defined. If /bnd = 0, 
the update is rejected. Otherwise, we randomly select a 
bond from the list and propose to remove it. The accep- 
tance probability for the update is 



P 



rmb 



IhndPAB 

{M+l)nB 



-/"(ra 



1 



(4.7) 



From Eqs. H4.6|l - H4.7() we realize that an optimal choice 
for Pab is based on the interaction potential between the 
bin centers 



An estimator for the bond contribution to the poten- 
tial energy, t/bonds, is obtained using a standard trick of 
replacing v{r) Xv{r) and then utilizing the identity 



1 dZ I 



(4.9) 



in accordance to which each Z-configuration is differen- 
tiated with respect to A, and then A is set to unity. The 
contribution to the derivative from the bond factors 14.4|l 
yields 



(4.10) 



-ui-RB^-RA). (4.8) 



where the sum is over all bonds in a current configuration. 



A. Worm updates within the diagrammatic bonds 

How does the presence of bonds affect the worm up- 
dates? The answer depends on the updating scenario. 
The easiest way is to work with the same updating pro- 
cedures we had without bonds. The only extra price is 
then in having simple constraints on the updates applica- 
bility to a given configuration. Namely, we require that 
all the beads being either deleted, or created, or shifted 
as a result of worm updates be free of bonds. Special 
care of bond factors has to be taken when interconvert- 
ing regular beads to worms [see Eqs. 14. 4|) and H4.5|l ]. 

The above requirement is not, in practice, as restrictive 
as one might think, since the probability of having no 
bonds at M consecutive beads is of order unity, with the 
proper choice of parameters discussed in the next Section. 

A brief qualitative discussion of the parameter Tc is in 
order here. Formally, can be as small as the size of 
the repulsive core of the potential, to guarantee the posi- 
tive definiteness of the second term in Eq. H4.3|l . It turns 
out, however, that this choice is not optimal, because 
in this case the total number of bonds per M consecu- 
tive beads of a world line may grow large (assuming that 
M is optimized in terms of the worm updates), and the 
above-mentioned condition of absence of bonds, in order 
for worm updates to be performed, may be satisfied very 
infrequently; consequently, the scheme may become inef- 
ficient. In the next Section we show that for helium, Tc 
is slightly larger than the radius of the first coordination 
shell. For best performance, both M and Tc should be 
simultaneously optimized. 

Finally, we note that it is possible to generalize the 
diagrammatic procedure to cases when the e"'^'^-' ex- 
ponential does not factor into pairwise terms. A great 
simplification here comes from the observation that for 
any particular choice of U{X), factorization does take 
place at least to the first approximation. The correcting 
terms are then of the form of close-to-unity three-bead, 
four-bead (and so forth) factors. The larger the number 
of beads in the correcting factor, the closer the correct- 
ing factor is to unity in terms of the powers of e. This 
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observation immediately suggests a recurrent (perturba- 
tive) scheme, which consists of ascribing correction fac- 
tors to the diagrams, depending on the number of beads 
in the bond-connected cluster. The description of such a 
scheme goes beyond the scope of the present paper, since 
we find that higher-order corrections can be safely ne- 
glected for the realistic choice of Vc discussed in the next 
Section. 



V. ADDITIONAL NOTES ON THE OPTIMAL 
ALGORITHM PARAMETERS FOR *HE. 

The separation radius for the diagrammatic expansion 
is determined from two conditions which ensure high per- 
formance of the algorithm. For efficient Swap updates 
one has to keep the length (in imaginary time) of the 
updated trajectory long enough to avoid proposing large 
displacements over short time periods. If ao is the inter- 
atomic distance, then the parameter M must satisfy the 
condition 



M 



(ao/2)'m 
2e 



which for "'He with oq = 3.5 A gives M ~ 25 assuming 
relatively small e = 5 x 10^'^ K~^. 

Since Swap and other updates are performed on trajec- 
tory pieces having no bonds on the corresponding time 
intervals, it is important to have also 



(^b„d(M)) 



1 



(5.2) 



where {lhnd{M)) is the average number of bonds on the 
trajectory interval of length Me. If this number is large, 
the simple updating strategy presented above will be- 
come inefficient due to small probability of fluctuations 
to the state with lhnd{M) = 0. In practice, the com- 
putational cost of calculating acceptance ratios in the 
strongly-correlated system is much larger than a simple 
check for the presence of bonds; thus, the above condition 
can be easily extended to {lhnd{M)) « 2. 

The number of bonds per time interval can be readily 
estimated from the pairwise interaction potential by us- 
ing a mean-field estimate for the chemical potential shift 

/>00 /"OO 

(Zbnd) « Men / d'^r v{r) g{r) « Men / d'^r v{r) , 



(5.3) 

where g(r) is the pair correlation function, which tends to 
1 in the r ^ oo limit. The last approximation in Eq. (|5.3|) 
is quite accurate for r greater than a few times cq. For 
the Aziz pair potential, Me — 0.125 i^^^, and particle 
number density n = 0.025 we have the condition 
(^bnd) ~ 2 satisfied for 



4.2 A 



(5.4) 



in three dimensions. This radius falls in between the first 
and second peaks of the pair correlation function, and. 



roughly speaking, includes the first coordination sphere 
of ^He atoms. Thus, on average, not more than 12 neigh- 
boring particles are contributing to the potential energy 
term in the exponent for the configuration weight of the 
trajectory. It should be stressed that the above set of 
control parameters is, of course, merely an approximate 
guideline. To give a concrete example, for simulations 
of liquid helium near the A-point at the saturated vapor 
pressure (SVP), we find that rc = 4 A and Me ~ 0.125 
is a reasonable choice. 

The advantage of using sophisticated forms of U{X), 
lies in the fact that one may achieve the same accuracy in 
evaluating the configuration weight, with smaller number 
of time slices (see Ref. for an exhaustive discussion of 
this aspect). Below we discuss the scheme suggested in 
Ref. 123, which takes into account the first derivatives 
of the interparticle potential (this is the scheme that we 
have adopted in all calculations presented below): 



UiX) 



j=2k 



- E 

j=2k+l 



4e 



9m 



(5.1) 

where 



(5.5) 



dv 



N.N 
»=1 jT^i 



N N.N 

nRj)-Efi^-E(E% 

i=l i=l ^ j^i 

(5.6) 

sums squares of forces acting on particles. The "force" 
term is very important in the repulsive region of the po- 
tential and close to the point where v{r) changes sign. In 
this region v{r) derivatives are large and their contribu- 
tion to the configuration weight may become comparable 
to the leading linear in e terms on some slices. 

Let us separate in Eq. H5.6|l contributions coming from 
forces between close (r^ < Tc) and distant (r^v > Tc) 
pairs (we label them with indexes 1 and 2 respectively), 

N . s2 N . s2 

^-E(Ef^O ^ ^ = 5:(Ef^0 ^ 



(5.7) 



If only Fi and F12 terms were present, the diagrammatic 
expansion would be easy to modify to include these terms 
into the consideration. The -Fi term is accounted for 
directly in the weight exponent by keeping records of 
short-range forces acting on particles and updating them 
accordingly. We assume that this procedure is always 
implemented and discuss below only how to deal with the 
force term in the diagrammatic expansion of the potential 
tail. 

The F12 term modifies the diagram weight and accep- 
tance ratios for the diagrammatic updates. Now, the rel- 
ative weight of configurations with and without a bond 
between distant (odd) beads a and b is given by the same 
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formula H4.4(l with 

... = f<r.,) + g(f«-fW).e- (5.8) 

Since the procedure of keeping track of forces acting on 
particles is standard for high-accuracy PIMC schemes, 
the required modifications of the scheme are minimal. 
[Note that short-range forces are present now in the ex- 
pression for the bond factor (|5.8|l . and thus their possible 
effects on bonds should be accounted for whenever these 
forces change]. 

Dealing with the F2 term is more cumbersome. Its ex- 
act treatment requires a solution of recursive relations for 
every chain of beads connected by the diagrammatic ex- 
pansion. At this point we notice that forces between dis- 
tant particles are orders of magnitude smaller than forces 
acting at short distances r < 2.6A and thus can be safely 
neglected. Indeed, for Tc = 4.5 A and e = 5 x 10^^ K^^, 
the F2 term can be estimated to be of order of 10~^ and 
thus has no measurable effect on the simulation results. 
After all, the interatomic potential is not even known 
with this accuracy. We conclude then that the diagram- 
matic expansion can be straightforwardly implemented 
for the high-accuracy scheme by omitting the F2 term in 
the configuration weight. 

Formally, due to extremely rare statistical fluctuations 
which bring two particles at very close distance and re- 
sult in large forces f acting on them, the sign of Uab in 
Eq. (|5.8() may change from positive to negative. This, in 
turn, will change the sign of the configuration weight if 
the corresponding bond is accepted. One can hardly clas- 
sify the possibility of such rare events as a "sign problem" 
because the average configuration sign will remain close 
to unity, and will not impair the algorithm efficiency. In 
practice, for Tc = 4 — 5 A and e < 0.01 K^^ the con- 
figuration sign simply never changes during the entire 
simulation. 

Finally, one may wonder if the inclusion of the F12 
term really helps to achieve better accuracy with smaller 
number of slices. If the desired accuracy does not ex- 
ceed three significant digits and Tc is kept larger than 
4 A then the answer is "No" . At this level of accu- 
racy, one may implement the diagrammatic expansion 
for the potential tail by ignoring the force term in the 
bond factor altogether, i.e. exactly as described in the 
previous section. Though we have implemented schemes 
with and without the F12 term, we did not detect any 
difference in final answers when using algorithm param- 
eters specified above. The Fi term was still kept in the 
exponent for more accurate evaluation of the short-range 
part. Alternatively, one may choose to deal with the 
short-range part using the pair-product approximation 
of Ceperley and PoUockf^ making use of tables for the 
two-particle density matrix calculated for an interaction 
potential which is identically zero at distances greater 
than Tc. 



VI. SIMULATION RESULTS 

In this Section, we present WA simulation results for 
bulk liquid "^He in two and three dimensions, with the 
aim of demonstrating that superfiuid properties and off- 
diagonal correlations can be calculated with the WA for 
very large system sizes, orders of magnitude larger than 
accessible to the conventional PIMC technology. We use 
the standard interatomic (Aziz) potential for helium, in 
an early form for consistency with other calculations.— 
The reason for our choice of illustrative system, is simply 
that the simulation of ''He in its condensed phase is a 
de facto test bench for new computational many-body 
techniques. 

As mentioned above, we have utilized the form for 
U{X) suggested in Ref. |23 for all the calculations for 
which results are presented here. In both two and three 
dimensions, we have observed convergence of the kinetic 
energy estimates (computed with the usual thermody- 
namic estimatoi"^) using a time step £=1/640 K~^; for all 
other quantities, four times a value of e can be used, and 
the estimates are seeing to coincide, within their statis- 
tical uncertainties, with those extrapolated to the e ^ 
limit. 



A. *He in two dimensions 

We start with the two dimensional case and extend 
the study of the superfluid-normal liquid (SF-N) transi- 
tion first carried out in Ref. 25, at a density n = 0.0432 
A^^; we consider system sizes up to hundred times larger 
than in the original study>2& In Fig. [S] we present data 
for the superfluid fraction (T) , for systems comprising 
N = 25, 200, and 2500 atoms. Since the SF-N transition 
in 2D is in the Kosterlitz-Thouless universality class^i 
with strong (logarithmic) finite-size corrections, a reliable 
extrapolation to the thermodynamic limit requires that 
simulations be performed for significantly different num- 
ber of particles. In the vicinity of the transition point, 
one may then employ the asymptotic (in the limit of large 
distances) vortex-pair renormalization group (RG) the- 
ory to fit the data. In our study, we used the same RG 
procedure as in RefiS^, which is based on the notion of 
the vortex core diameter, d (as a short-range cut-off for 
RG equations), and vortex energy, Ec, to control vortex 
density at distance d. Only data in the narrow vicinity 
of the transition point 0.65 K <T < 0.8 K were used in 
the fitting procedure. Our estimates for the values of the 
fitting parameters are (i = 8.8±0.5Afor the vortex core 
diameter, and i?=2.18 ± 0.04 K for the vortex energy, 
which lead to an estimate for the critical temperature 
rc=0.653±0.010 K. The vortex diameter (roughly twice 
the interatomic distance) turns out to be comparable to 
the linear system size L for N = 25, clearly showing that 
the use of the asymptotic RG analysis is questionable for 
such a small system size. This fact also explains why our 
result for is significantly different from the previous 
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FIG. 5: Superfluid fraction ps{T) computed for 2D *He on 
systems with different numbers A'^ of ''He atoms. The system 
density is n = 0.0432 A~^. Dashed lines represent fits to 
the numerical data (in the critical region) obtained using the 
procedure illustrated in Rof.^''. The leftmost dashed line is 
the extrapolation to the infinite system. Open squares show 
results obtained in Refi— for the same system, with N—25. 




r(A) 

FIG. 6: (Color online) . One-particle density matrix computed 
for 2D *He at a density n — 0.0432 for a system of 

200 atoms, at r=0.675 K (upper curve) and r=1.0 K (lower 
curve). Statistical errors on the curves are very small, and 
not shown for clarity. In the inset, we present data (on a log- 
log scale) for the N=2500 system at T=0.675 K, with clear 
signatures of the Kosterlitz-Thouless behavior in the vicinity 
of the critical point. [Reproduced from Ref. ITU] 



estimate, 0.72±0.02 K, deduced from the = 25 data.^ 
In Fig. El we show resuhs for the single-particle density 
matrix n{r). For 2D hehum, this quantity is expected to 
decay to zero at all finite temperatures. In the normal 
phase, and far from the critical point, the decay is expo- 
nential. At the critical point, the decay is described by a 
slow power- law n{r) - l/r^/^. The same law should be 
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FIG. 7: (Color online). One-particle density matrix n(r) close 
to the SVP critical point at T — 2.14: K for two system 
sizes A'^ = 64 (filled squares) and A'' — 2048 (open squares). 
The solid line is the theoretical prediction based on long- 
wavelength phase fluctuations, Eq. 16.111 . 



observed in the vicinity of the critical point up to expo- 
nentially large distances. At low T, the exponent in the 
power law approaches zero. The data in Fig.|^are in line 
with these expectations; quantitatively, the slope of the 
n(r) curve on the log-log plot shown in the inset is close 
to one quarter. 

We wish to conclude this subsection by giving an ex- 
ample of typical (without any extensive optimization) al- 
gorithm parameters used for the two-dimensional helium. 
For the T = 1 K, fi = -1 K, and A^ = 25 system with 
P = 200, M = 40, Co = 7.5, and = 4.05 A the mea- 
sured acceptance probabilities were Pop = 0.47, Pel = 
0.43, Pin = 0.13, P™ = 0.44, Pad = 0.43, P,, = 0.59, 
P,w = 0.33, P„b = 0.58 P^b = 0.42. 



B. 



He in three dimensions 



As mentioned in the Introduction, simulations of bulk 
3D liquid helium were among the first remarkably suc- 
cessful applications of the PIMC methodii However, pre- 
vious predictions made for the superfluid properties were 
never at the same level of accuracy as for energetic or 
structural properties, for reasons mentioned in the Intro- 
duction. In this subsection, we show how the WA, based 
on updates described above, eliminates the shortcomings 
of the existing PIMC method, by allowing simulations of 
several thousand atoms with sufflcient precision to deter- 
mine, for example, the critical temperature of the SF-N 
transition at the saturated vapor pressure (SVP) with 
accuracy of three significant digits. 

In Fig.UI we show our data for the density matrix n(r) 
in the vicinity of the critical point, at the saturated vapor 
pressure (at a temperature T=2.14 K and at a density 
n=0.02198 A-3) for system sizes A^ = 64 and A^ = 2048. 
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Though the data for small and large system sizes appear 
nearly identical in the main plot, the crucial difference is 
clearly seen in the inset. The iV = 64 curve terminates 
right after the first coordination shell oscillation; the best 
estimate that can be obtained of the condensate fraction 
Uo (namely, the asymptotic value to which n{r) should 
plateau at long distances) from this set of data alone, 
would be about 0.045. Obviously, the same coordination 
shell oscillation prohibits a fortiori any reliable finite- 
size scaling for smaller system sizes, e.g. using series 
N = 16, 32, 64. 

In contrast, the N = 2048 system is large enough to see 
the effect of long-wavelength hydrodynamic phase fluc- 
tuations. The Bogoliubov expression for the asymptotic 
behavior of n{r) at large distances in the superfluid, is 
given by 



n{r) = Ho exp 



T 



SirXripsr 



(6.1) 



Since ps is calculated independently, the shape of the 
density matrix decay is fixed. The condensate fraction 
controls only the overall normalization of the theoretical 
curve, and this allows precise extrapolation of the data 
to the thermodynamic limit. An example of such ex- 
trapolation is shown in Fig. [7| by the solid line, which 
predicts uq = 0.024(1) for the condensate fraction Uo at 
T = 2.14 K — nearly a factor of two smaller than the 
0.045 estimate obtained on a 64-atom system. 

The hydrodynamic correction to the tail of n(r) is less 
important at low temperature and for large values of ps . 
For comparison, in Fig. |S1 we present data for a system 
of = 1024 atoms, at a temperature T — 1 K {in this 
case, the density is n=0. 02184 A"'^). In this temperature 
range, smaller system sizes can be used to obtain reliable 
thermodynamic estimates of Hq. Our estimate for tiq at 
T=l K is 0.082±0.002. This is consistent with the ex- 
isting PIMC estimate (0.07±0.01 at T=1.2 K, from Ref. 
0), obtained on a system of 64 ''He atoms, but some- 
what above the most recent r=0 estimate (0.069±0.005, 
Ref. I28j) . Although the difference is only slightly greater 
than the combined statistical uncertainties, it should be 
noted that finite temperature calculations are unbiased, 
whereas ground state calculations are based on an input 
trial wave function. While it is in principle possible to 
remove the variational bias associated to the trial wave 
function, this may be difficult a goal to achieve in prac- 
tice. 

Having access to system sizes which allow asymptotic 
hydrodynamic description is a necessary condition for 
determining critical parameters using finite-size scaling 
techniques (see, e.g., Ref l2^. The idea is to consider 
quantities which are determined by system properties at 
the largest scales, becoming scale-invariant at the critical 
point. For small deviations from criticality, (5 — > 0, the 
dependence on system size for such quantities is given by 
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FIG. 8: One-particle density matrix n(r) for A'^ = 1024 par- 
ticles at T = 1 K and SVP pressure. The solid line is the 
theoretical prediction based on the long-wave phase fluctua- 
tions, Eq. 16. H . 



where /fl(x) and gn^x) are the corresponding universal 
scaling functions {gii{x) is analytic at a; = 0) and ^(S) is 
the correlation length which diverges at the critical point 
as ^ oc 5^'^. For the U(l) universality class in 3D, the 
best numerical estimate currently available for the corre- 
lation length exponent is = 0.6717.^" The intersection 
of R{L, 6) curves for different system sizes provides very 
accurate and unbiased estimates of critical parameters. 

Previous attempts to determine Tc from the scale- 
invariance of R{L) = p^L (Josephson relation) have 
failed Z Though the authors of Ref. 7 correctly argue that 
"The statistical uncertainty in the data is too large to ac- 
curately determine Tc from the crossing of these [scaling] 
two curves" , it seems that coordination shell oscillations 
also contribute to several intersections in the 1.6 K < T < 
2.4 K interval. 

In Fig. O we show the temperature dependence of the 
superfluid fraction for various system sizes. After extrap- 
olation to the thermodynamic limit, there is nearly per- 
fect agreement between the numerical and experimental 
results. In order to determine the transition temperature 
Tc, we perform flnite-size scaling analysis of 



R(L,T) = 



2\npsL _ 



T 



R[L,5)^ fn{L/a5))^gn{5L^/'^) 



(6.2) 



where 'W={Wx,Wy,Wz) is the winding number The 
raw data are shown in Fig 1101 As an independent 
check, we also draw a horizontal line at the known U(l)- 
universality class value for winding number fluctuations 
at the critical point From the intersection of scal- 
ing curves (which is seen to take place at the univer- 
sal value, within the statistical uncertainties), we find 
Tc — 2.193(6). The difference between this prediction 
and the experimental value, 2.177 K, is very small, i.e., 
simulations of the superfluid density dependence on sys- 
tem size in the vicinity of the critical point, do allow us 
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FIG. 9: (Color online). Superfluid fraction ps(T) as a function 
of temperature at SVP, computed for different system sizes, 
namely A'^ = 64 (filled circles), = 128 (open circles), A'^ = 
256 (filled squares), iV = 512 (diamonds), = 1024 (triangles 
down), and A' = 2048 (triangles left). The solid line is the 
experimental curve. 
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FIG. 10: (Color online). Finite-size scaling plot for 
2XnpsL/T = {W^)/3 at SVP. The solid line is the U(l) uni- 
versality class value of 0.516(1). 
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FIG. 11: (Color online). Condensate fraction at the saturated 
vapor pressure. The dashed line is used to guide the eye. 
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FIG. 12: (Color online). Liquid density as a function of chem- 
ical potential at low temperature T — 0.25 K. The slope of 
the solid line is deduced from the simulated values of com- 
pressibility averaged over four points shown in the plot. The 
critical value of the chemical potential deduced from this plot 
is estimated as /lo = 0.06 ± 0.04. 




to calculate Tc with a relative accuracy of better than 
0.5%. There is no reason to expect the Aziz pair po- 
tential to reproduce Tc with much better accuracy, as 
it was not optimized for this purpose. It is only a very 
good approximation to the true interatomic potential for 
*He which, in reality, includes irreducible forces acting 
between three and more particlesj^^ 

Condensate density has been consistently one the most 
difficult properties to compute and measure for helium. 
Shown in Fig. ^] are available experimental datai^ along 
with the previous PIMC results? and new estimates ob- 
tained in the present study, extrapolated to the thermo- 
dynamic limit as explained above. With new technology 



we substantially reduce theoretical uncertainties on pre- 
dictions of Uo at the SVP. 

We conclude by mentioning that we have also obtained 
estimates for all other standard thermodynamic quanti- 
ties, such as the kinetic energy, pair correlation function 
etc. Our data are generally consistent with those of exist- 
ing calculations; specifically, our T=l K results are indis- 
tinguishable, within statistical uncertainties, from those 
yielded b y nu merically exact ground state methods (see, 
e.g., Ref.jH. 
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VII. DISCUSSION AND CONCLUSIONS 

We have introduced a new Worm algorithm, afford- 
ing an accurate PIMC study of strongly correlated Bose 
systems. Illustrative results of numerical simulations of 
the superfluid transition in liquid ^He have been pre- 
sented, for system sizes two orders of magnitude larger 
than what is accessible to conventional PIMC. It should 
be stressed that such an advance cannot be simply at- 
tributed to the availability of faster computing facilities 
than back in the days when the first PIMC simulations of 
liquid ^He were carried out. Rather, the WA decisively 
overcomes the most important limitation of conventional 
PIMC, namely the exponential inefficiency with which 
long permutation cycles are sampled, in the thermody- 
namic limit (a limitation acknowledged by practitioners 
of PIMC*^). 

We have also described a procedure, based on ideas of 
Diagrammatic Monte Carlo, which allows one to enhance 
significantly the scalability of the computational scheme, 
without compromising on the accuracy of the calculation. 

The new methodology has already been applied to the 
study of the supersolid phase of heliumfi^ for which ac- 
cess to large system sizes is crucial; it can also be ex- 
pected to have broad impact on a wide variety of strongly 
correlated quantum many-body systems. It should be 
mentioned, that the efficiency with which long permu- 
tation cycles can be sampled using the WA, also signifi- 
cantly impacts the convergence of calculations of the su- 
perfluid properties of finite systems, such as quantum 
clusters 

There are other advantages to this new method, chiefly 
the fact that it is fully grand canonical, and that allows 
for the calculation of the Matsubara Green function, a 
quantity that can not e computed with any other existing 
QMC technique (in continuous space). An immediate 
application of these last two aspects, is the calculation 
of chemical potentials and excitation gaps. Consider, for 
instance, the calculation of the chemical potential at the 
liquid-solid transition line at low temperature, /io. The 
study of G{k — 0, t) in the solid phase can be found in 
Ref.lM 

The value of fio can be determined from the curve 
and known freezing density for the liquid which is at 
n = 0.02599 A~^. In Fig. [T^lwe show data points for 
four values of the chemical potential and the position of 
the freezing density. All simulations were performed at 
T = 0.25 K, and with the particle number N « 800. 
Apart from average density we also calculated statistics 
of particle number fluctuations to obtain the compress- 
ibility of the system from k = dn/dn = {{N-{N)y)/TV. 
It provides an independent check for consistency and con- 
vergence of the data. In Fig.^jthe slope of the solid line 



is obtained from the average value of n for all data points 
shown. From this set of data we deduce that the critical 
value of the chemical potential as 

Aio = 0.06 ± 0.04 . (7.1) 
l.OIr 1 

1^ ^ 
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FIG. 13: Zero-momentum Green function for the superfluid 
state of N ^ 1000 atoms at T = 1 K, p = -7.35 K, and 
density n = 0.02184 A~"^ . It is normalized to unity at the 
origin. Note the small vertical scale: The r-dependence is a 
finite-size effect. 

Though one can deduce this parameter from the ground 
state energy as a function of particle number using fi = 
dE/dN « E{N +1) ~ E{N), such simulations can not 
be performed for large system sizes with accuracy signif- 
icantly better then few K. 

The zero-momentum Matsubara Green function G{k — 
0, r) for the superfluid state at T = 1 is shown in Fig. 1131 
In the macroscopic limit, G(fc = 0, r) is r- independent, 
being equal to the total number of condensate particles. 
We show this plot to support our previous claim that in 
the superfluid state the two worm ends perform a random 
walk on large distances, and come close within the reach 
of the Close and Remove updates once per sweep. 
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